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In cosmic ray air showers, the muon lateral separation from the center of the shower is a measure 
of the transverse momentum that the muon parent acquired in the cosmic ray interaction. IceCube 
has observed cosmic ray interactions that produce muons laterally separated by up to 400 m from the 
shower core, a factor of 6 larger distance than previous measurements. These muons originate in high 
Pt (2 - 15 GeV/c) interactions from the incident cosmic ray, or high-energy secondary interactions. 
The best fit to the separation distribution includes a power law component, indicating the presence 
of a hard pr component that can be described by perturbative quantum chromodynamics. However, 
the rates and the zenith angle distributions of these events are not well reproduced with the cosmic 
ray models tested here, even those that include charm interactions. This discrepancy may be 
explained by a larger fraction of kaons and charmed particles than is currently incorporated in the 
simulations. 



I. INTRODUCTION 

There have been many attempts to measure the cosmic 
ray composition at energies around and above the knee 
of the spectrum (~ PeV) [l], 0. At these energies, di- 
rect measurements by balloon and satellite experiments 
have very limited statistics. Ground based experiments 
rely on indirect measurements using observables such as 
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the ratio of the measured electromagnetic energy to the 
number of muons These analyses are dependent on 

phenomenological calculations and simulations to relate 
the muon observations to an inferred composition; the re- 
sult can be sensitive to the assumed hadronic interaction 
models @. 

Studies of high-energy (> 1 TeV) muons with under- 
ground detectors have been an important part of this 
effort. The muons are produced early in air showers and 
probe the initial shower development Q. Two classes of 
muons are generally considered. "Conventional" muons 
come from pion and kaon decays, while "prompt" muons 
come from the decays of particles containing heavy 
quarks, mostly charm. Conventional muons dominate at 
TeV energies, but, at energies above 100 TeV, prompt 
muons are expected to dominate Q. The resulting 
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change in the slope of the muon energy spectrum has 
not yet been observed ||. 

Studies of isolated muons, far from the shower core, can 
help understand the uncertainties due to phenomenolog- 
ical models. Muon separations greater than about 30 m 
are largely due to the transverse momentum, px, im- 
parted to the muon by its parent. For px > 2 GeV/c, 
these interactions can be described in the context of per- 
turbative quantum chromodynamics (pQCD). Data from 
RHIC, the Tevatron, and the LHC are in quite good 
agreement with modern fixed order plus next-to-leading 
log calculations @. These experimental studies give us 
some confidence in pQCD calculations for air showers. 

Experimentally, the transition from soft interactions 
(i.e. those with px < 2 GeV/c that are not describable 
in pQCD) to hard interactions is visible as a transition 
from a pt spectrum that falls off exponentially to one 
that follows a power law. At low pt the spectrum follows 
exp (—px/T), with T « 220 MeV/c for pions (somewhat 
higher for kaons and protons) [lOj. At higher px, the 
spectrum falls as 1/(1 + px/po) n , where one fit found 
n=13. 0to° 5 and p = 1.9±g;f GeV/c |Il|. The transition 
is around 2 GeV/c for pions. This spectral change should 
be visible in the lateral separation distribution. 

The MACRO detector has previously measured the lat- 
eral separation between muons in air showers for pri mary 
energies ranging roughly from 10 4 GeV to 10 6 GeV [12j. 
Buried under 3.8 km water equivalent of rock, MACRO 
has a minimum muon energy of about 1.3 TeV. MACRO 
measured muon pair separations out to a distance of 
about 65 m. Their simulations verified the linear rela- 
tionship between px and separation (with a small offset 
due to multiple scattering of the muons) out to a px of 
1.2 GeV/c, below the expected transition to the pQCD 
regime. 

The muon px is related to the separation of the muon 
from the shower core by 



dx ~ 



pxHc 
cos{0) 



(1) 



where dx is the perpendicular separation between the 
muon and the shower core, H is the interaction height of 
the primary, H/ cos(#) is the path length of the shower 
to the ground at a zenith angle 0, and is the energy 
of the muon at generation. The interaction height of the 
parent of the muon is assumed to be synonymous with the 
primary interaction height, and the energy of the muon 
at generation is well approximated by its energy at the 
surface of the Earth. 

In addition to the initial px, muons can separate from 
the shower core when they bend in the Earth's magnetic 
field or multiple scatter in the ice above the detector. 
However, for the muon energies and separations consid- 
ered here the gyroradius is on the order of 20,000 km and 
multiple scattering is similarly negligible; the initial px 
is the dominant effect producing the separation. 

A selection of muons with large transverse separation is 
biased toward events produced at high altitudes, high px, 



and low energy. The detector geometry imposes a min- 
imum energy; most of the muons will naturally be near 
this threshold. In the competition between altitude and 
Px, the atmospheric density decreases exponentially with 
increasing altitude, while the px should fall more slowly, 
only as a power law. Although widely separated muons 
are biased toward high-altitude interactions, they still 
allow for studies of transverse momentum. Events with 
large zenith angles are also preferred, since the muons 
are given more time to separate. Of course, all of these 
factors should be appropriately modeled in the Monte 
Carlo. 

In the IceCube neutrino telescope [IH, muons are de- 
tected with a 1 km 3 array of optical sensors buried in 
the Antarctic ice at depths between 1450 and 2450 m; 
individual muons studied in this analysis must have an 
energy at the Earth's surface of at least 400 GeV to reach 
the detector. The 125 m horizontal spacing between Ice- 
Cube strings serves as a rough threshold for the minimum 
resolvable separation. For vertical muons with an inter- 
action height of 50 km and an energy of 1 TeV, this corre- 
sponds to a px of 2.5 GeV/c. The interaction height and 
muon energy vary from shower to shower, so the event- 
by-event uncertainty in px approaches a factor of 2 if we 
assume average values for both. 

The muon energy and cosmic ray interaction height 
for air showers that produce muons in the detector de- 
pend on the zenith angle. Figure [T] shows the interaction 
height as a function of true zenith angle for all showers 
at sea level for DPMJET simulation (see Sees. IIVI and 
IVIII for a full description of the simulation). This de- 
pendence arises in part because the column depth of the 
atmosphere as well as the primary energy for air showers 
increases with zenith angle, leading to higher interaction 
heights for horizontal showers. The interaction height for 
showers with a muon at least 135 m from the shower core 
is also shown; these showers interact much higher in the 
atmosphere. Figure [2] shows a fit to both the minimum 
and average energy of muons at the surface of the earth 
as a function of zenith angle, as well as the minimum 
energy calculated assuming continuous energy loss along 
the track (dE/dx). The energies shown are for the muon 
with the largest perpendicular distance from the shower 
core. 

The zenith angle has an impact on shower develop- 
ment. The 1450 m of ice above the detector shield it from 
vertical muons with energies less than about 400 GeV. 
For inclined showers, several effects come into play. The 
distance between the target and detector rises, giving the 
muon more time to separate from the shower core. How- 
ever, the slant depth also increases, raising the muon en- 
ergy threshold roughly exponentially with the ice thick- 
ness. For average dE/dx energy loss, the minimum en- 
ergy at the Earth's surface is given by fl4| 
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where D is the depth of the detector and a and b are 
constants that describe the energy loss of a muon in ice 
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FIG. 1. (Color online). The interaction height for all DPM- 
JET simulated showers. Distributions for simulated showers 
with a true maximum muon separation less than 135 m and 
greater than 135 m are also shown. 
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where the three terms are for muons from pions, kaons, 
and charmed particles, respectively. Charmed hadrons 
decay very quickly, leading to a flatter distribution in 
zenith angle. This angular difference can be used to sep- 
arate prompt muons from conventional muons; Fig. [3] 
shows the zenith angle of cosmic ray muons with ener- 
gies of ~2 TeV produced by pion, kaon, or charm inter- 
actions. The fraction of muons from charm interactions 
increases for high zenith angles and higher energies. 
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FIG. 3. (Color online). The cosine of the zenith angle of DP- 
MJET simulated muons produced by pion, kaon, or charmed 
particles interactions. The curves have been normalized to 
the peak bin and the muons energy is 2 TeV (±10%). 



FIG. 2. (Color online). The true energy at the surface of 
the Earth for the muon furthest from the shower core for 
simulated Sibyll shower events that pass all selection criteria 
versus zenith angle. Also shown are fits of the minimum and 
average energies. The values calculated using Eq. [2] are also 
shown. 



(with values of 0.177 GeV/mwe and 0.209 x 10~ 3 /mwe, 
respectively) (l4| . 

The zenith angle distribution also depends on the par- 
ents of the muons. At TeV energies, most muons origi- 
nate from pions and kaons that decay before they inter- 
act. The probability of decay increases for larger zenith 
angles, because the pions and kaons spend more of their 
livetime at higher altitudes where the target density is 
lower so they are less likely to interact. The muon flux 

is u m 



This paper extends the MACRO muon lateral separa- 
tion measurements out to a separation of 400 m, well into 
the pQCD regime, using 335 days of data collected with 
the partially completed IceCube detector. The following 
sections give a description of the IceCube detector and 
an overview of the analysis. The simulation is described 
in section IIVI and the background reduction is discussed 
in section [V] The resulting distributions are discussed in 
sections [Vl] and I VIII 



II. THE ICECUBE DETECTOR 

IceCube is a 1 km 3 underground neutrino telescope lo- 
cated at the South Pole. 5,160 digital optical modules 
(DOMs) on 86 vertical strings detect Cherenkov radia- 
tion from charged particles traversing the Antarctic ice 
[H|. The DOMs are located between 1450 and 2450 m 
below the ice surface. Each DOM consists of a 25 cm 
photomultiplier tube [16j plus associated digitization and 
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calibration electronics [17[ , all in a 35 cm diameter pres- 
sure vessel. 

The primary energy and pseudorapidity of the muon 
with the maximum separation in simulated IceCube cos- 
mic ray events that survive all selection criteria for this 
analysis are shown in Fig. |4] The majority of muons 
studied here are produced in the far-foward region. 
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FIG. 4. Top: The true primary energy of Sibyll simulated 
shower events that pass all selection criteria. The center of 
mass energy is shown on the top axis. Bottom: The pseudora- 
pidity of the muon with the maximum transverse separation 
in Sibyll simulated shower events that pass all selection crite- 
ria. 

This analysis uses data taken from May 20, 2009 to 
May 30, 2010 when the array was partially complete, 
with 59 of the 86 strings deployed. Each DOM triggers 
internally at a threshold that gives it an 85% efficiency 
for single photoelectrons; data is sent to the surface if the 
trigger occurs in coincidence with the DOM's nearest- 
neighbor or next-to-nearest-neighbor. There, a surface 
trigger selects events where 8 DOMs triggered within 
10 fj,s. The detector triggered at a rate of about 1600 Hz; 
the bulk of these events were from cosmic ray muons. 



III. ANALYSIS OVERVIEW 

This analysis is sensitive to a shower with a bundle 
of low pt muons that are contained within a few 10s 
of meters of the shower core, plus an isolated high pt 
muon separated at least 125 m from the shower core. The 
bundle gives a good reconstruction of the core location 
and primary zenith angle. 

Each cosmic ray event is reconstructed with a two- 
track hypothesis. One track reconstructs the center of 
the muon bundle, and the other track reconstructs the 
muon with the largest lateral separation, which we call 
the laterally separated muon (LS muon). The perpen- 
dicular separation between the two tracks at the point of 
closest approach to the center of IceCube is defined as 
the lateral separation of muons in a cosmic ray shower 



(d,T in Eq. [T]). For showers without a distinct separation 
between the bundle core and the LS muon, this mea- 
surement can be thought of as the lateral extent of the 
cosmic ray air shower. Figure [5] shows the simulated true 
LS muon distance from the shower core at the initial filter 
level and after applying all selection criteria. For show- 
ers with separations greater than ^100 m, the spectral 
behavior of the iron and proton showers are very similar. 




Maximum Muon Separation From Shower Core [m] 

FIG. 5. (Color online). The Sibyll simulated true separation 
between the LS muon and the shower core for all cosmic ray 
showers, showers with proton primaries, and iron primaries at 
the EHE filter level as well as the distribution after applying 
all selection criteria (described in Sec. \V\ . 



The background for this measurement is cosmic ray 
air showers without LS muons. There are two distinct 
types of background. The first is single showers from 
cosmic rays that appear as a single track in the detec- 
tor and can be eliminated by requiring two high quality 
reconstructed tracks. Multiple independent coincident 
showers ("double showers") constitute the second type 
of background. The IceCube 59-string configuration is 
large enough that the rate of simultaneous cosmic rays 
events is significant. Muon bundles from two (or more) 
uncorrelated air showers can deposit light in the detec- 
tor within the event window, producing two (or multiple) 
separated tracks. The rate of double showers scales with 
the single shower rate and falls exponentially with the 
number of coincident showers, making double cosmic ray 
events the only significant background. The number of 
double showers can be reduced by requiring that the two 
reconstructed tracks arrive at roughly the same time (i.e. 
within ±450 ns of each other) and from the same direc- 
tion. However, an irreducible background remains from 
showers that arrive simulataneously from the same direc- 
tion. This number is very small compared to the signal 
and can be measured by studying the data in the "off- 
time" window. Off-time events have tracks that arrive 
between 450 ns and 1350 ns of each other, after applying 
all other selection criteria. 
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IV. SHOWER SIMULATION 

Cosmic ray air showers are generated using the COR- 
SIKA simulation program. The simulated cosmic ray 
spectrum was based on the Horandel polygonato model 
[ill with energies between 600 GeV and 10 11 GeV. The 
showers are generated using the Sibyll 2.1 [2(J, DPM- 
JET 2.55 [lH, and QGSJETOlc [H| hadronic interaction 
models that are used to develop the shower in the atmo- 
sphere. QGSJET and DPMJET include charm interac- 
tions, but the current version of Sibyll does not (see 23 1 
for a detailed comparison of these models). The muons 
were propagated through the ice by the Muon Monte 
Carlo simulation package [24| . The propagation of light 
through the ice and the response of the electronics, as 
well as the trigger, were simulated with IceCube software. 
All simulations were done with the CORSIKA MSIS-90-E 
atmospheres paramaterized for South Pole atmospheres 
measured in March, July, October, and December of 1997 

Double shower events are simulated by combining sin- 
gle shower events with a random separation in time. The 
detector response to the double showers is simulated as 
one event. Events are counted as double showers even 
if they do not produce enough hits in the detector to 
reconstruct a track. This approach is conservative and 
results in an overestimation of the double shower rate 
for this analysis. This is acceptable because the double 
shower rate is measured using the data and the simula- 
tion serves only as a guide for event identification and 
development of selection criteria, and to calculate event 
rates at early cut levels before the two track reconstruc- 
tion is performed. 

Signal and background distributions are drawn from 
the same dataset, where showers with muons more than 
100m from the shower core that reach IceCube depths are 
considered LS muon signal. Roughly 23 days of simulated 
livetime have been generated with Sibyll, 11 days with 
QGSJET, and 7 days with DPMJET. 

The majority of the figures shown in this paper are 
made with Sibyll because it has the largest simulated live- 
time. For most of the figures there is no significant differ- 
ence between the different hadronic interaction models, 
so the lines for the other models are omitted for clair- 
ity. Distributions that required detailed knowledge of 
muon production inside the air shower were drawn from 
dedicated simulation that recorded the full interaction 
history of all muons. Due to computational limitations 
only the two hadronic models that included charm in- 
teractions, QGSJET and DPMJET, were used for the 
dedicated simulation. 

A test simulation with maximally varied ice properties 
showed a 20% increase in simulated event rate. However, 
this uncertainty is smaller than the variation in event rate 
for different hadronic interaction models, so this system- 
atic is not included. Generating air showers is compu- 
tationally intensive and computational resources are the 
primary limit on the livetime of generated simulation. 



Nonetheless, for this analysis the statistics generated are 
more than sufficient for comparison to data. 



V. ANALYSIS 

This analysis uses the extremely high energy (EHE) fil- 
ter stream where events are required to generate at least 
630 photoelectrons j2~6| ; the filter output rate is about 
1.4 Hz. The photoelectron requirement corresponds to a 
primary energy threshold of about 1000 GeV, nearly an 
order of magnitude below the minimum energy of events 
that survive all the selection criteria (Fig. [4]) and is suffi- 
ciently low to avoid bias in the final selection. The EHE 
stream is sensitive to events from all directions and al- 
lows comparison of event rates as a function of zenith 
angle. After the EHE filtering, 6.4 x 10 7 events were left 
in the data sample. Additional reconstructions were run 
on the remaining events. 



A. Reconstructing LS Muon Events 

LS muon events have a unique dual topology: a bundle 
of low pt muons that make up the core of the shower, 
and a laterally separated muon with the same direction 
and timing. These two components are reconstructed 
separately. 

An initial, fast reconstruction based on a linear rela- 
tionship between the arrival times and light wavcfront 
gives the approximate direction and location of the bun- 
dle track. The hits are rotated into a plane perpendicular 
to the first guess track. The rotated hits are sorted into 
two sets using the k-means clustering algorithm, an al- 
gorithm that sorts the hits into two clusters according to 
their closeness to the mean of the cluster [glj • The larger 
set of hits is assumed to be from the muon bundle and is 
reconstructed with a maximum-likelihood function that 
accounts for the arrival time of the Cherenkov photons 
and the scattering of light in the ice [28| . 

Hits that belong to the LS muon can be identified by 
their timing relative to the reconstructed bundle track. 
Since the LS muon arrives in the detector at roughly the 
same time as the bundle but at least a hundred meters 
from the shower core, its light has a much earlier arrival 
time than light could propagation from the shower core. 
Hits are considered LS muon hits if their arrival time is 
more than 100 ns earlier than the expectation for light 
from the bundle track. Additionally, LS muon hits are 
also required to be more than 90 m from the bundle fit to 
reduce the contamination of hits belonging to the bundle. 
The values for the timing and separation were chosen 
to minimize miscategorization of hits and maximize the 
number of events where two tracks can successfully be 
reconstructed. 

Next, the process is iterated to increase the accuracy of 
the reconstruction. First, the hits are resorted according 
to their perpendicular distance from the reconstructed 
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LS muon track. Any hits that are more than 100 m from 
the LS muon track are used to reconstruct a new bun- 
dle track. LS muon hits are selected based on the same 
timing and separation values used previously but rela- 
tive to the new bundle track. For most events, the new 
reconstructed tracks are not very different, but iterating 
helps eliminate the few cases where the initial clustering 
algorithm did not perform well. 

Figure [B] shows the cumulative fraction of events as 
a function of the space angle between the reconstructed 
and true track for Sibyll simulation. The reconstruction 
algorithms are able to reconstruct the direction of the 
bundle to within 4.0° and the LS muon to within 5.6° 
for 68% of the events. Figure [7] shows the resolution of 
the measured separation between the two tracks for Sibyll 
simulation. The muon with the largest separation from 
the bundle center is used for the true value of the LS 
muon. The separation is measured with a resolution of 
~30 m. 




5 10 15 20 

Space Angle between Generated and Reconstructed Track 



FIG. 6. (Color online). The angular resolution of the bundle 
and LS muon reconstructions for Sibyll simulation at the final 
selection level. 
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FIG. 7. The separation resolution of the bundle and LS 
muon reconstructions for Sibyll simulation at the final selec- 
tion level. 



B. Background Reduction 

Background reduction was studied with Sibyll since it 
has the largest livetime. Table fl] shows the number of 
events passing each selection level. As mentioned in Sec- 
tion IIV| the number of double showers in the simula- 
tion is overestimated. This is used to calculate the total 
background event rates at early selection levels before 
the two-track reconstruction was performed. Once the 
two track reconstruction is available, rates are shown for 
double showers calculated using off-time events in data. 

After initial filtering, only events where the bundle and 
LS muon track reconstructions both succeeded are re- 
tained. This reduced the background from single show- 
ers by roughly two orders of magnitude while retaining 
more than 300,000 of the signal LS muon events. The 
main reason signal events fail to reconstruct is there are 
not enough hits in the LS muon track; the average num- 
ber of hit DOMs from the LS muon is 20 (out of -3,500 
DOMs on 59 strings). The number of hit DOMs in the 
bundle and LS muon groups are shown in Fig. M 




Numbered hit rjiOMs inLS n 



FIG. 8. (Color online). The number of hit DOMs in the bun- 
dle group (top) and the LS muon group (bottom) for data and 
Sibyll simulation at the final selection level. The simulation 
curve is normalized to the data. 



1. Double Showers 

The background from double showers is reduced by re- 
quiring that the direction of the two reconstructed tracks 
agree to within 5° and that they arrive at the point of 
closest approach to the detector center within ±450 ns 
of each other. An irreducible background remains from 
double showers that happen to come from the same direc- 
tion at the same time. Requiring the two tracks originate 
from the same direction and arrive at the same time re- 
duced the rate of double showers by a factor of 1500, 
while reducing the simulated signal by a factor of two. 
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TABLE I. Number of events passing each selection level in 335 days for data, the sum of all background (EBG), LS muons 
(dr > 100), single showers and double showers estimated from Sibyll simulation (sim), and double showers estimated from 
off-time data. The uncertainties shown are statistical. 



Selection Criterion 


Data 


EBG 






LS Muons 




Single 






Double (sim) 


Double 


EHE Filter 


63,649,372 


5.47 ± 0.02 


X 


10' 


1.77 ± 0.03 


X 


10 b 


4.31 ± 0.02 


X 


10' 


9.79 ± 0.08 x 


10 b 




Reco. Successful 


3,804,388 


2.87 ± 0.01 


X 


10 6 


3.24 ± 0.02 


X 


10 s 


3.94 ± 0.02 


X 


10 5 


2.15 ± 0.01 x 


10 6 




|AT| < 450 ns 


723,592 


6.99 ± 0.03 


X 


10 s 


3.12 ± 0.02 


X 


10 5 


2.66 ± 0.02 


X 


10 5 






120,384 


A$ < 5° 


327,962 


3.24 ± 0.02 


X 


10 5 


1.84 ± 0.02 


X 


10 5 


1.38 ± 0.01 


X 


10 5 






1,445 


Single Likelihood > 7.5 


69,877 


7.54 ± 0.10 


X 


10 4 


7.02 ± 0.10 


X 


10 4 


4.14 ± 0.24 


X 


10 3 






1,020 


d T > 135 m 


56,463 


5.94 ± 0.09 


X 


10 4 


5.74 ± 0.10 


X 


10 4 


1.08 ± 0.13 


X 


10 3 






905 


LS Muon nDOMs > 8 


38,966 


4.71 ± 0.08 


X 


10 4 


4.57 ± 0.08 


X 


10 4 


6.35 ± 0.96 


X 


10 2 






695 


LS Muon nStrings > 2 


34,754 


4.52 ± 0.08 


X 


10 4 


4.42 ± 0.08 


X 


10 4 


5.62 ± 0.90 


X 


io- 






456 



2. Single Showers 



3. Angular Resolution 



Elimination of single showers relies on the topological 
difference between single and double tracks. Single tracks 
are well-reconstructed by the likelihood functions made 
for a single track hypothesis while events with LS muons 
are not. Figure [9] shows the log of the output of the 
likelihood function for a single track reconstruction, done 
using all the hits (bundle and LS muon) , for single shower 
background, LS muon signal, and data. Events that are 
good fits to the single track reconstruction hypothesis 
have a lower value on this plot. Requiring a likelihood 
greater than 7.5 reduced the background by a factor of 
30 and retains 40% of the signal. 




Data 
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FIG. 9. (Color online). The log likelihood for reconstruction 
with a single track hypothesis for data, all simulated Sibyll 
showers, LS muons, single showers, and double showers esti- 
mated from off-time data. 



Next, events with a dr less than 135 m are removed. 
The closer the two tracks are, the more difficult it is for 
the reconstruction algorithm to separate them. Remov- 
ing these closer tracks decreased the single shower back- 
ground by a factor of four, while only removing 20% of 
signal events; significantly improving the signal to back- 
ground ratio. 



After applying selection criteria to remove single and 
double showers, the predicted number of background 
events was about l/30th the number of LS muon events. 
This was more than sufficient for measuring the LS muon 
separation spectrum, but some of the LS muon events 
still had poorly reconstructed directions. These LS muon 
tracks mostly triggered DOMs on a single string, making 
it difficult to reconstruct azimuth information. These 
events were eliminated by requiring that the LS muon 
hits occur on at least 3 IceCube strings. Also, the LS 
muon track was required to be robust, triggering more 
than 8 DOMs. These cuts select high quality events im- 
proving the angular resolution for the LS muon track 
from 7.3° to 5.6° while only removing 20% of the signal 
events. 



VI. RESULTS 

After applying all the selection criteria, 34,754 events 
remain in data. Of those, the expected number of ran- 
dom double showers (based on off-time data rates) is 
456. The number of predicted events depends on the 
interaction model. The Sibyll simulation predicts 45,200 
events, 98% of which are LS muon events, while simu- 
lations with QGSJET and DPMJET predict 57,300 and 
27,500 events, respectively. 

Figure [TU] shows the lateral distribution of data events 
that pass all selection criteria. The expectation from sim- 
ulated LS muons is also shown and is in good agreement 
with the data. Additionally, the separation of double 
showers for off-time data events is shown. The separa- 
tion distribution of double showers is flat, consistent with 
random coincidences (modulo detector edge effects). Fig- 
ure QT] shows the ratio of the simulation to the data as 
a function of lateral distribution in the region where LS 
muon signal is dominant. A linear fit to the ratios did not 
find a statistically significant difference in the slopes for 
the three simulation distributions, indicating they are all 
a reasonably good match for the data. Figure shows 
the lateral distribution of muons in cosmic ray showers 
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FIG. 10. (Color online). The separation between the LS 
muon and bundle track after applying all selection criteria 
for data, simulation with the Sibyll, DPMJET, and QGSJET 
interaction models, as well as double showers estimated from 
off-time data. 



FIG. 12. (Color online). The separation between the LS muon 
and bundle track for data events after applying all selection 
criteria and subtracting the off-time events for three zenith 
bands. 



x 2 / ndf 45 / 21 




Separation between Bundle and LS Muon [m] Separation between Bundle and LS Muon [m] 



FIG. 11. (Color online). The ratio of simulation to data 
versus the separation between the LS muon and bundle track. 

for three different zenith slices after subtracting the sep- 
aration distribution from the double showers. The shape 
of the separation distribution does not show any strong 
dependence on zenith angle, indicating the separation is 
not a strong effect of propagation length. 

A. Distribution At Sea Level 

The separation distribution shown in Fig. 1101 includes 
the effects of detector efficiency. To remove these, we 
applied a bin-by-bin efficiency correction calculated from 
the ratio of the simulated separation distributions at sea 
level and after applying all selection criteria for Sibyll, 
QGSJET, and DPMJET combined. The resulting distri- 
bution is shown in Fig. 1131 

If the larger separations are due to large transverse 
momenta, the expected distribution would follow an ex- 



FIG. 13. (Color online). The LS muon data distribution at 
sea level, along with the best fit parameters for the function 
described in the text. The exponential part of the fit is plotted 
as a dotted red line and the power law is shown as a dashed 
blue line. 



ponential that transitions to a power law at large sep- 
arations. To test this, the separation distribution was 
fit with a function that is the sum of exp(A + Bx) and 
10 c (l +x/400) n . A, B, C, and n were allowed to vary. 
The resulting fit, shown in Fig. [T31 has a best fit with a 
transition at 308 m to a power law that falls off as sep- 
aration to the -11. This composite fit has a x 2 /DOF of 
45/21, compared to a x 2 /DOF of 68.6/23 for a purely 
exponential function. 

For comparison, we have applied the same fit to the 
Monte Carlo models. For the sum of Sibyll, QGSJET, 
and DPMJET, the x 2 / D OF was 19.1/21 for the two com- 
ponent fit, compared to a x 2 /DOF of 958/23 for an ex- 
ponential fit. 
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B. Zenith Angle Distribution 

A comparison of the arrival angle shows significant dis- 
agreement between simulation and data. Figure [TJ] shows 
the cosine of the reconstructed zenith angle of the bundle 
for events that survive all selection criteria. Sibyll and 
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FIG. 14. (Color online). The cosine of the reconstructed bun- 
dle track after applying all selection criteria for data, simu- 
lated showers using the Sibyll, QGSJET, and DPMJET in- 
teraction models, as well as double showers estimated from 
off-time data (scaled by 10). 

QGSJET overpredict the event rate at high zenith angles 
and underpredict the rate for the more vertical events, 
while DPMJET is a better match to the data at least 
at high zenith angles. There is no significant difference 
between the QGSJET and Sibyll distributions. Figure 
[T5l shows the ratio of simulation to data as a function 
of zenith angle. A linear fit to the ratio showed that 
DPMJET had the lowest slope, indicating it is the best 
match of the three simulations. A Kolmogorov-Smirnov 
(KS) test of the simulated distributions against the data 
finds that none of the simulations are a good match. The 
highest correlation is found for the DPMJET distribu- 
tion, but the probability of it being drawn from the same 
distribution as the data was only 5 x 10~ 12 . 



VII. DISCUSSION 

Figure [16] shows the minimum px for a separation of 
135 m calculated with the minimum and average muon 
energies (shown in Fig. [5]) and the interaction height for 
showers with distant muons (from Fig. [1]). 

DPMJET is a better match for the data at high zenith 
angles where the px is expected to be slightly higher. 
However, even the DPMJET simulations underpredict 
the number of events at low zenith angles that also have 
high pt values. 

The different approaches to calculating px taken by 
the various interaction models may contribute to the dis- 
crepancy in zenith angle. Sibyll 2.1 generates particle 
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FIG. 15. (Color online). The ratio of simulation to data ver- 
sus the cosine of the zenith angle of the reconstructed bundle 
track after applying all selection criteria. 
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FIG. 16. (Color online). The minimum muon transverse mo- 
mentum of simulated shower events that pass all selection 
criteria for different energy parameterizations as a function of 
zenith angle. The interaction height comes from Fig. [T] 



Pt using two approaches [20j. In soft (low px) interac- 
tions, quark-antiquark pairs are given a px that follows 
a Gaussian distribution, with a mean of 0.3, 0.45, or 0.6 
+ 0.08 log 10 ^/s/30 GeV GeV, for light quarks, strange 
quarks, and diquarks respectively; diquarks lead to the 
production of baryons. Hard interactions are simulated 
using lowest order pQCD, with the GRV structure func- 
tions, including a saturation correction. The minimum 
px to qualify as a hard interaction increases as a func- 
tion of energy. Charmed particles are not included in this 
version of Sibyll. 

DPMJET 2.55 also simulates both hard and soft inter- 
actions, with the hard interactions based on lowest order 
pQCD, using simple, phenomenological structure func- 
tions [23]. Charmed particles are produced both by soft 
processes, either at the ends of soft sea chains or inside 
the chain decay, and by "hard processes at the ends of 
hard and semi hard chains (minijets)" [3(|. The latter 
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mechanism is modeled using pQCD, and the soft pro- 
cesses use phenomenological models. All possible charm 
particles are simulated. 

QGSJETOlc simulates interactions using cross sections 
calculated with an eikonal that is the sum of contribu- 
tions from soft and hard processes 22|. Charmed par- 



ticles are produced when a charm quark-antiquark pair 
from the vacuum is coupled to the hadronizing strings of 
soft or semi-hard jets. Only the lightest charmed meson 
and baryon are generated (the D and A c ). A comparison 
of predicted cross section and px distributions with ex- 
perimental measurements showed reasonable agreement 
at low pt [ll| . 

While all of these simulations include a hard compo- 
nent, it is interesting to note that DPMJET is the only 
model that includes a hard component of charmed par- 
ticles, and it is also the model that agrees best with the 
data. 

None of the current Monte Carlo codes include bottom 
quark production. At these energies, bottom quark pro- 
duction is only a small contributor to the overall muon 
flux. However, at RHIC, bottom quark production is a 
significant contributor to the prompt lepton flux at high 
Pt, contributing more than 30% of the prompt lepton 
production with px > 2 GeV/c in proton- proton colli- 
sions at a center of mass energy of 200 GeV [32| . Similar 
fractions are seen in 7 TeV proton-proton collisions at 
the LHC 33[. This will increase the expected prompt 
high pt muon flux by about 50%. 

The difference in zenith angle distribution may also 
indicate a larger fraction of muons from heavier parents. 
According to Eq. [31 the muon distribution becomes flat- 
ter as the parents become heavier; this behavior can be 
seen in Fig. [3] The flatter data distribution is consistent 
with a larger fraction of muons from kaons or charmed 
particles than is predicted in the simulation. A dedi- 
cated simulation that recorded the parents of all muons 
was performed for QGSJET and DPMJET. A full sim- 
ulation of detector response was not possible, so events 
similar to events in the final sample were selected by re- 
quiring the muon energy to be greater than the minimum 
muon energy shown in Fig. [2] and that the event have a 
muon at least 100 m from the shower core. Figure 1171 
shows the fraction of muons with maximum separation 
from the bundle produced by pions, kaons, and charm 
particles in the shower. DPMJET has a higher frac- 
tion of muons produced by kaons at the zenith angles 
where the agreement with data is best. This suggests 
that muon production by kaons plays an important role 
in these events. Experimental measurements have shown 
a tendency for the kaon fraction to increase with target 
mass; various estimates of the kaon-pion ratio vary by as 
much as 20% [H. 

The discrepancy in zenith angle may also be related 
to differences in the composition of the cosmic rays. Fig- 
ures fTS]- [2UJ show the ratios of the simulated cosine of the 
bundle zenith angle to data for Sibyll, QGSJET, and DP- 
MJET for several different primary compositions. Sibyll 
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FIG. 17. (Color online). The fractions of muons with max- 
imum separation from the shower core from pions, kaons, 
and charm for DPMJET (closed shapes) and QGSJET (open 
shapes) for a sample of events similar to the final events in 
the data (see text for details). 



and QGSJET show a flatter ratio for protons than for 
heavier elements. A KS test shows that the proton-only 
distributions are the best match to the data distribu- 
tion for these two simulations. Contrastingly, the ratios 
with DPMJET simulations are all relatively flat with rel- 
atively similar KS probabilities. 

This suggests an interplay between kaon and charm 
abundance and composition. It is easier for proton pri- 
maries to produce high px muons because all of their en- 
ergy is concentrated in one particle. Sibyll and QGSJET, 
which have fewer kaons and (no) charmed particles, can 
only reproduce the data distributions with the cosmic 
rays most likely to produce high pt muons: protons. 
However, DPMJET simulation, with its greater fraction 
of kaons and charmed particles, is closer to the data dis- 
tributions for all cosmic ray compositions. 
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FIG. 18. (Color online). The ratio of Sibyll simulation and 
data as a function of cosine of the reconstructed bundle zenith 
for different cosmic ray primary compositions. 
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FIG. 19. (Color online). The ratio of QGSJET simulation 
and data as a function of cosine of the reconstructed bundle 
zenith for different cosmic ray primary compositions. 
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FIG. 20. (Color online). The ratio of DPMJET simulation 
and data as a function of cosine of the reconstructed bundle 
zenith for different cosmic ray primary compositions. 



The better agreement with data for protons simulated 
with Sibyll and QGSJET may indicate a dependence on 
the relative abundances of cosmic ray primaries. A cos- 
mic ray composition model with lighter composition be- 
tween 1 to 10 PeV could also improve the agreement in 
zenith angle for these models. 

The zenith angle distribution may also depend on the 
inelasticity in ttN and KN interactions. If pions and 
kaons lose less energy than predicted when they inter- 
act, then the high-energy muon flux may be enhanced, 
especially near the vertical direction. 



VIII. CONCLUSION 

IceCube has observed 35,000 muons with lateral sep- 
arations between 135 m and 400 m; this corresponds to 
a transverse momentum in the range of 2 to 15 GeV/c. 
The separation distribution is poorly fit by an exponen- 
tial distribution with a x 2 / D0F of 68.6/23. The fit im- 
proves when a power law component is included to a 
X 2 /DOF of 45/21, as expected from pQCD. However, the 
zenith angle distribution of the muons is unexpectedly 
flat, even when including the decay of charmed particles, 
and is poorly modeled by current simulations. This may 
be caused by an underproduction of kaons and charmed 
particles in the simulation. 

IceCube has demonstrated the capability to resolve lat- 
erally separated muons in air showers. Future analyses 
with modified and improved simulations could lead to an 
estimate of muon parent ratios as well as a measurement 
of the transverse momentum spectrum in cosmic ray air 
showers. 
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